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^»,i , Abstract. Coalescing compact object binaries consisting of black holes and/or 

■ Neutron stars are a prime target for ground-based gravitational wave detectors. 

jH ' This article reviews the status of numerical simulations of these systems, with an 

emphasis on recent progress. 
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1. Introduction 

Inspiraling and coalescing compact object binaries consisting of black holes and/or 
Neutron stars are of high importance to understand gravity and matter in extreme 
conditions. The strong gravitational fields require full nonlinear general relativity to 
describe them, and so compact object binaries allow exploration of curved space-time 
in the nonlinear and highly dynamic regime. The matter density inside Neutron stars 
is above nuclear density. Neutron stars therefore provide an avenue to study matter 
at super-nuclear densities, and compact object mergers allow such studies even in 
dynamic situations when the two objects collide. 

Compact object binaries are therefore at the centre of attention for several 
scientific disciplines: The study of black hole binaries (BH-BH) sheds light on 
properties of general relativity in the genuinely nonlinear, dynamic regime (in a 
clean vacuum environment), ranging from black hole kicks (e.g. [1]) to the topological 
structure of event horizons [2-4] . Simulations of binaries involving one or two Neutron 
stars (BH-NS, NS-NS) [5, 6] elucidate the connection between these systems and short 
gamma ray bursts. Finally, gravitational waves emitted by compact object binaries are 
the most promising source for gravitational wave detectors like Advanced LIGO [7- 
9], Advanced Virgo [10,11] and KAGRA [12]. Gravitational wave detectors require 
accurate waveform models and information about electromagnetic counterparts to 
optimally search for gravitational waves via matched filtering ("event detection"), 
and to extract the maximum amount of information about the source of gravitational 
waves, once such waves have been detected ("parameter estimation"). 

Direct numerical simulations of the full Einstein equations are an important 
cornerstone in the study of coalescing compact object binaries. These simulations 
have continued to progress swiftly, resulting in a large number of notable recent 
results. This article gives an overview of the current status of numerical simulations of 
compact object binaries. Our particular focus lies on recent advances, which have not 
yet been incorporated into longer or more specialized review articles (like [6, 13, 14]) 
or books (like [15]). Numerical simulations of higher-dimensional gravity (e.g. [16, 17]) 
and alternative theories of gravity (e.g. [18,19]) are advancing rapidly, but because 
of length limitations we shall restrict our attention to the "classical" compact object 
binaries (BH-BH, BH-NS, NS-NS) in four space-time dimensions in standard general 
relativity, the scenario of most direct relevance to gravitational wave detectors. 

This article is organized as follows: Section 2 reviews BH-BH simulations, 
starting with numerical methods, recent advances, applications to gravitational wave 
science, and ending with BH-BH simulations embedded in gaseous or electromagnetic 
environments. Section 3 deals with binaries with Neutron stars (BH-NS, NS-NS). We 
conclude in Sec. 4 with some thoughts about the future of the field. 

2. Black Hole — Black Hole Binaries 

2.1. Numerical methods 

Since the numerical relativity breakthroughs in 2005 [20-22], several frameworks 
have emerged to simulate inspiraling and merging BH-BH binaries. The major 
differentiating factors are the choices for the formulation of Einstein's equations (either 
the BSSNOK system [23-25] or the Generalized Harmonic (GH) system [20,26-28]) 
and the choice of the numerical evolution algorithm, finite-differences (FD) or pscudo- 
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Table 1. Ingredients into a BH-BH simulation, and the most common choices 
for the two major frameworks in use to solve BH-BH spacetimes. 



spectral methods. Pretorius' first simulations employed the combination GH+FD [20, 
27], and this code has since been applied to a variety of physical scenarios including 
high energy collisions of BH-BH [29] , eccentric BH-BH [30] and studies of gravity in 
different space-time dimensions [17,31]. 

However, in regard to the main focus of this section, BH-BH simulations of 
direct relevance to gravitational wave detectors, the combinations BSSNOK-I-FD and 
GH-|-Spectral have emerged as the leading frameworks, and so we shall describe these 
in some detail. The BSSNOK+FD framework is buih around the BSSNOK [23-25] 
evolution equations coupled with finite-difference (FD) approximation schemes. The 
GH-|-Spectral framework utilizes the generalized harmonic (GH) equations [20,26- 
28] and multi-domain spectral methods. These frameworks differ in a large number 
of the individual elements that are required to compute a gravitational waveform 
with numerical relativity. Table 1 lists the main choices that are made within the 
two frameworks for setting initial data, performing the evolution and analyzing the 
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resulting data. Besides formulation and numerical algorithm used for the evolution 
equations, these choices include formulation and numerical methods for the initial 
data, diagnostics tools, and factors that determine the physics of the scenario under 
consideration (e.g. low orbital eccentricity). While Table 1 conveys the broad picture, 
it must necessarily simplify and omit details. For instance, Pretorius' breakthrough 
work [20] does not fit into these categories^. 

The choices shown in Table 1 are interrelated and only certain combinations 
are feasible. In particular, one cannot simply exchange specific elements between 
BSSNOK+FD and GH+Spectral. Let us illustrate some of these dependencies: 
Moving puncture evolutions performed in the BSSNOK+FD framework require very 
specific gauge conditions (Gamma-driver and 1+log shcing [45-47]) to keep the 
punctures stable during the evolution. Because the BH's are not excised, moving 
puncture evolutions require initial data covering the entire initial-data hypersurface. 
Therefore, a moving puncture evolution cannot be started from the excision initial- 
data used for the spectral evolutions§. 

Puncture data covers the entire initial-data hypersurface, so it seems that 
GH-|-Spectral evolutions could use puncture data. However, this is also non-trivial: 
Singularity excision requires careful control over the location and shape of the excision 
boundary [49, 56]. This is more difficult when the black hole horizons change on short 
time-scales, in particular during merger|| and early in an evolution when the coordinate 
location of the black holes and the gauge settles down. Conformal thin sandwich initial 
data (as used in the GH+Spectral approach) reduces significantly such transients and 
is therefore more suitable for GH+Spectral than puncture initial data. 

Trumpet initial data reduces initial transient apparent horizon motion [57] 
compared to standard puncture initial data, however, so far trumpet initial data not 
been used regularly for production BH-BH simulations. 

To summarize some key properties of the two BH~BH frameworks: 

• BSSNOK+FD simulations are generally considered more robust and "easier" 
because of the lack of black hole excision and the larger expertise of how 
to stabilize finite-difference methods. Moreover, code-infrastructure [58-60], 
including adaptivc-mesh-refinement (Carpet [61,62]) and apparent horizon 
finders [63], are publicly available. Several independent research groups use these 
tools and have obtained a tremendous wealth of results (e.g. black hole kicks). 
The major codes are listed in Table 1 and more details are available in the related 
proceeding contributions describing the NINJA collaboration [64]. 

• Spectral methods, in contrast, require BH excision and arc very sensitive to the 
presence of ill-defined elements in the problem (e.g. outer boundary conditions, 
handling of inter-domain boundaries where grids of different resolution touch, 
or filtering). Significant care and expertise is needed to analyze and control all 
relevant aspects of the evolution, and only one spectral BH-BH code exists: The 
Spectral Einstein Code SpEC [65], developed by the SXS collaboration between 
Cornell, Caltech, CITA and Washington State University. 

t Pretorius used the generalized harmonic evolution system with an adaptive-mesh-refinement finite- 
difference code, starting from initial data that contained nonsingular balls of scalar field that collapsed 
to BH's early in the evolution. 

§ Proposals to fill-in the BH interiors of excision initial-data [53-55] have not been systematically 
pursued. 

II Control of the excision boundary was one of the major obstacles in obtaining BH-BH mergers with 
GH+Spectral, cf. [49,56]. 
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• The key advantage of spectral methods is their fast (exponential) convergence. 
Therefore, high accuracy can be obtained with comparatively low computational 
cost. The high accuracy/lower cost permeates the entire "simulation pipeline": 
Lower numerical noise (due to higher accuracy) makes it easier to perform certain 
fits that are required in iterative eccentricity removal (see Sec. 2.2.2); the low 
computational cost makes it possible to perform longer simulations; the clean 
waveforms extracted at finite radius are more amenable to extrapolation to obtain 
the asymptotic waveform at infinitely large distance; high accuracy also makes 
it easier to compute and analyze quantities that depend on higher derivatives of 
the evolved fields, for instance black hole vortices and tendices [66,67]. 

• To date, many more simulations have been performed with BSSNOK+FD than 
with GH+Spectral, but the GH+Spectral simulations arc longer and more 
accurate. Both techniques have found their place, depending on the precise 
needs of the science to be studied. However, this conventional wisdom slowly 
changes, as the BSSN+FD simulations become more accurate and achieve longer 
simulations, whereas the GH+Spectral simulations become more robust and the 
number of GH+Spectral simulations increases [64]. 

While simulations work successfully, we caution the reader that both frameworks 
utilize hcuristically determined ingredients: Gauge conditions and constraint damping 
terms contain parameters chosen by trial and error, and mesh structures are 
tuned based on user experience. Therefore, parameter choices working in one 
region of parameter space may require adjustments for other regions (e.g. [68-70]). 
Furthermore, constraint damping parameters impact the accuracy and quality of 
the simulation beyond their main purpose of preserving the constraints [71]. The 
presence of user-tuned parameters has further implications: On the negative side, it is 
not guaranteed that the current techniques work in unexplored regions of parameter 
space. On the positive side, it is conceivable -even likely- that further tuning of code 
parameters will enhance accuracy and efficiency of future simulations compared to 
todays runs. 

2.2. Recent advances 

Given the rather mature state of BH-BH simulations, recent advances are incremental, 
pushing and refining capabilities, and confirming assumptions that were made in 
previous years to achieve fast progress. 

2.2.1. New records Let us start with some "new records": With increasing mass- 
ratio q = Mi/M2 > 1, fully numerical simulations become more challenging because 
the computational cost increases (the small black hole must be resolved) , and because 
experience of code tuning from comparable mass binaries is no longer applicable. Two 
groups pushed the mass-ratio of BH-BH simulations up to g= 100: Sperhake et al [72] 
examine head-on collisions of a small black hole with a large black hole. Lousto and 
Zlochower [73] perform a numerical simulation of the last two orbits of a BH-BH 
binary of mass-ratio q= 100. Both these calculations are very impressive, requiring 
not only a large amount of CPU resources, but also improvements and tuning of the 
computational infrastructure (e.g. location of mesh-refinement regions). The RIT 
group [74-76] extracts BH-trajectories from numerical simulations of mass-ratio up 
to g = 100. These trajectories are then used in perturbative calculations to compute 
the waveforms of high mass-ratio binaries. 
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Turning to spin, the long-standing bound 

]|2=X<0.93 (1) 

was recently broken. Certain black hole properties vary significantly between x = 0.93 
and X = 1 and therefore angular momentum may not be the most useful measure 
of extremality. For instance, the rotational energy of a black hole with x = 0.93 is 
only 59% of the rotational energy of a maximally rotating black hole. When taking 
rotational energy as the metric, the bound (1) is far from extremality. 

Equation (1) represents the maximal spin that is achievable with puncture initial 
data [77], and therefore, all BH-BH simulations within the BSSN-OK framework 
are limited by this bound. BH spins higher than the bound (1) are achieved with 
the more general initial-data formalism utilized within the GH-|-Spectral framework. 
Lovelace et al [78] construct initial data and perform short evolutions of BH-BH 
binaries with spins of 0.97. More recently, Lovelace et al [56, 79] compute complete 
inspiral/merger/ringdown simulations for equal mass BH-BH binaries with aligned 
spins of 0.97 and anti-aligned spins of 0.95. These simulations require delicate control 
of the excision boundaries and demonstrate that the SpEC code has significantly 
matured. 

The X = 0.97 aligned spin simulation [56] also carries another new record: It is 
presently the longest published complete BH-BH simulation, lasting about 25 orbits 
before merger and ringdown, for an evolution time of 7000M. The longest incomplete 
simulation (i.e. inspiral only) so far is presented by Le Tiec et al. [80], lasting 34 
orbits and 11,000M. All these extremely long simulations were obtained with SpEC. 
The simulations presented in Ref. [80] range up to mass-ratio 8, all in excess of 20 
orbits, and are used to compare periastron advance in BH-BH binaries with several 
analytical approximation schemes. 

The last year has also witnessed another increase in the largest known black hole 
kick from inspiralling BH-BH^. To remind the readers, in 2007 the "BH super-kick 
configuration" was discovered [82-84], with BH kicks close to 4000km/sec. In 2011, 
Lousto and Zlochower [1] demonstrated that tilting the BH spins from the super- 
kick configuration toward partial alignment with the orbital angular momentum can 
increase the BH kicks to almost 5000km/sec. Partial alignment of the BH spins with 
the orbital angular momentum allows the BH's to spiral closer to each other, where 
the higher velocities enhance the BH kick. 

2.2.2. Improved techniques Eccentricity removal has been extended to precessing 
binaries: Isolated BH-BH binaries originating from binary stars are expected to have 
vanishing orbital eccentricity when they enter LIGO's frequency band [85,86]. To 
achieve low eccentricity in a numerical simulation requires careful choices for orbital 
frequency fig and radial velocity Vr of the individual black holes (or, equivalently, 
tangential and radial linear momentum of the BH's). These parameters can be 
chosen based on post-Newtonian information, either in closed form or by integrating 
post-Newtonian ordinary differential equations for two point-masses starting at large 
separation. The post-Newtonian coordinates and velocities are then used in the 
construction of BH-BH initial-data. This approach [42, 87] is computationally 
inexpensive and achieves eccentricities of a few 0.001 for low-spin and comparable 

1 Interactions of high- velocity black holes result in even larger kicks [81] . 
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mass BH-BH, and somewhat higher eccentricity for high spins and unequal mass BH- 
BH. 

One can also adjust orbital parameters iteratively [43,88]: One begins with a 
reasonable first guess for flo and Vr (perhaps based on post- Newtonian information), 
evolves for about two orbits, analyzes the orbital trajectories, and then adjusts Hq 
and Vr- With the adjusted initial orbital parameters, one constructs a new BH-BH 
initial data set, and performs a new evolution. This technique is computationally more 
expensive but achieves eccentricities as small as ~ 10^^. Recently, Buonnano et al [44] 
extended this technique to precessing BH-BH systems. Ref. [44] performs a post- 
Newtonian analysis demonstrating that spin impacts iterative eccentricity removal 
only for eccentricities 10""' . . . 10"'^ (the exact bound depends on the BH-spins and 
the initial black hole separation). Furthermore, Ref. [44] demonstrates that e ^ 10~'* 
can indeed be reached for several BH-BH configurations with dimensionless spins of 
0.5. This work uses the GH-|-Spectral approach with the SpEC code. Tichy et al [89] 
proposed an alternative iterative technique within the BSSNOK-I-FD framework and 
reach eccentricities of a few xlO~^. Very recently Piirrer et al. [90] study iterative 
eccentricity removal for moving puncture simulations based on the gravitational 
waveforms. 

Cauchy Characteristic extraction (CCE) [91-94] extracts certain data from 
a standard 3-1-1 evolution code (as described in Sec. 2.1) and evolves these data 
with a separate characteristic code to future null infinity. At future null infinity, 
gravitational radiation is unambiguously defined and so this technique promises gauge 
invariant waveforms, improving on the more widely used technique to extrapolate 
finite-radius waveforms to infinite extraction radius [95]. In recent papers, Babiuc et 
al [93,94] improve the PITT CCE code, and perform careful convergence tests. This 
code is publicly available as part of the Einstein Toolkit [60,96]. With the recent 
improvements and public availability of the CCE code, I expect rapid increase in the 
use of this post-processing tool. While CCE is gauge- invariant, it requires initialization 
at the begin of the numerical simulation. Bishop et al [97] investigate the impact of 
different initialization strategies and found some impact onto the CCE waveform. 
Furthermore, CCE cannot remove errors that were introduced in the underlying 
3-|-l numerical evolution. For instance, if the outer boundary conditions of the 3-1-1 
evolution admit unphysical incoming radiation (or refiected outgoing radiation), then 
such radiation will appear in the CCE waveform. 

2.2.3. Code validation and consistency checks Several papers confirmed that BH-BH 
simulations indeed work as expected thus further validating the numerical techniques: 

Owen [98,99] (extending work by Campanelli et al [100]) analyzes carefully 
a SpEC simulation of an equal-mass, non-spinning BH-BH. He confirmed with a 
multipolar analysis that the remnant black hole settles down to Kerr with the correct 
quasi-normal mode falloffs. The combination of gauge-invariant quantitites with the 
high accuracy of the SpEC code allows Owen to confirm agreement to many significant 
digits, usually to better than 1 part in 10^ and sometimes to better than 1 part in 
10^. Refs. [98, 100] also confirm that the space-time of a black hole merger approaches 
Petrov type D at late times after the BH-BH merger. 

Hinder et al [101] perform a careful analysis of the asymptotic fall-off of the 
Newman Penrose scalars, revisiting [102]. The fall-off rates agree with expectations, 

*n - 4r. (2) 
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where n = 0, . . . 4 labels the individual Ncwman-Penrosc scalars, vl'o, . . . , ^4. 

The NINJA-2 project [103] collected about 40 numerical waveforms, with 
some configuration contributed multiple times by different codes. These duplicate 
waveforms allow consistency checks between waveforms computed independently by 
different codes. Overlap calculations between (2,2) modes of hybridized waveforms 
show disagreements broadly in line with the expected errors of the different numerical 
codes and hybridization procedures. Comparisons of higher order modes have not 
yet been performed and are planned for future work. For details, see the separate 
NINJA-2 contribution to the Amaldi proceedings [64]. 

2.2.4- Future improvements Several recent directions of research may have an impact 
on accuracy and efficiency of BH-BH simulations in the future: The conformal 
Z4-system [104] combines features from BSSNOK with the constraint damping of 
the generalized harmonic equations, resulting in improved suppression of constraint 
violations (compared to BSSNOK). Witek et al [105] developed a generalized BSSNOK 
system, which encompases BSSNOK as a special case. They find that their extension 
improves numerical behavior of BSSN. Existing BSSNOK codes can easily be adopted 
to either the conformal Z4 system or the generalized BSSNOK system. Bona et al [106] 
presented a Lagrangian for the Z4 system, which facilitates the use of symplectic 
integrators for Einstein's equations. 

Improved efficiency or reduced wall-clock run-time are in dire need, given that a 
high-quality BH-BH simulation requires months to complete. At least two approaches 
are under development: Several groups [96, 107, 108] work on porting numerical 
relativity codes to graphical-processing- units (CPUs), usually within the CUDA 
framework, opening up the possibility that in the future CPUs may accelerate BH-BH 
simulations by a significant factor. Lau et al [109,110] explore novel time-stepping 
algorithms which circumvent the Courant timcstcp limit, and promise the ability to 
take significantly larger timesteps than current codes (see also [111]). 

2. 3. BH-BH waveforms for gravitational wave astronomy 

The primary motivation for the intense activity in BH-BH research lies in gravitational 
wave astronomy. Gravitational wave detectors require reasonably accurate waveform 
templates to detect GWs, and more accurate waveforms to extract detailed knowledge 
about the source properties (parameter estimation). As discussed in detail in the 
review [112], one first performs BH-BH simulations at discrete points in parameter 
space. One then constructs analytical waveform models (continuous in parameters) 
from these BH-BH simulations, which are used for GW data-analysis purposes. This 
sequence of steps has been carried out several times [113-120], based on different 
numerical simulations, different types of analytical waveform models, and for different 
regions of parameter space (again, we refer to Ref. [112] for details). 

The NINJA projects [64, 103, 121] give a good sense of the rate of progress of 
BH-BH simulations for GW-modeling. The initial NINJA project [121] in 2008 
consisted of about 20 numerical relativity simulations lasting on average about 12 GW 
cycles before merger, without any length- and accuracy-requirements. The current 
NINJA-2 project [64, 103] collected about 40 waveforms with on average about 20 
GW cycles. NINJA-2 also expands the coverage of parameter space by including more 
spinning waveforms, the simulations arc more accurate, and one carefully attaches 
post-Newtonian inspirals to the numerical waveforms. The second major ongoing 
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effort, the NR-AR collaboration [122] is assembling a yet larger number of waveforms 
of average length of about 30 GW cycles, and with yet higher accuracy than NINJA-2. 
Each one of these efforts takes about two years to complete, demonstrating the high 
complexity of computing, validating, and collecting high-quality BH-BH waveforms. 

One obstacle toward computation of high-quality waveforms in a shortage of 
experienced researchers capable of running BH-BH codes and capable of improving 
the efficiency, robustness and automation of the codes. Furthermore, there were 
unexpected difficulties in scaling BSSN-OK simulations to greater length and higher 
accuracy and for the SpEC code to obtain robust and automatic mergers. Finally, 
computational resources are also constrained: As a rule of thumb, for "easy" parameter 
choices (moderate spins ^ 0.5, moderate mass-ratios ^ 2, moderate length, ^ 10 
orbits), a single state-of-the-art BH-BH simulation requires on the order of 100,000 
CPU-hours+. This CPU-cost is of course dependent on the length T/M of the 
simulation (measured in units of the total mass A/). For given symmetric mass- 
ratio ;/, initial orbital frequency Vli or number N of orbits to merger, lowest order 
post-Newtonian expressions [123] yield: 

Zl;« A^-i(M^,)-8/3 (3) 

M 256 ^ ' ' ^ ' 

|^«5^3/^(27^iV)«/^ (4) 

Halving the initial orbital frequency fli or doubling the number of orbits N increases 
T by a factor ^ 6 or '^ 3, respectively. The increase in CPU-cost is even higher, to 
preserve phase-accuracy over the longer inspiral. Higher mass-ratio and higher spins 
increase the CPU-cost further. 

Equations (3) and (4) basically force a trade-off between the length of each 
simulation and the number of simulations that can be performed with limited 
CPU resources. Reasonable accuracy for event-detection can be achieved with 
~ 10 orbits [124,125], however, optimal parameter extraction requires numerical 
simulations starting at much lower initial frequency [126-128]. The low starting 
frequency is necessary because the accuracy of the hybrid waveform is primarily limited 
by the errors of the 3.5-th order post-Newtonian waveforms that are attached before 
the start of the numerical simulation. 

2.3.1. Exploration of parameter space The entire parameter space for BH-BH 
binaries is nine-dimensional: mass-ratio, two spin- vectors, and two parameters related 
to eccentricity (eccentricity and phase at periastron) . Essentially all efforts to explore 
this parameter space so far have focused on non-eccentric binaries, and even for non- 
eccentric binaries, only the following low-dimensional subspaces have been covered in 
detail (see the NINJA-2 contribution [64] for details and references): 

• Non-spinning binaries with mass-ratio 1 < 9 < 10. 

• Equal mass binaries with equal spins parallel to the orbital angular momentum. 

• Circular, non-precessing binaries form a three dimensional parameter space (mass- 
ratio and spin-magnitudes of the two spins parallel with the orbital angular 
momentum). This space has not been covered extensively yet, with most efforts 
having been focused on the two one-dimensional subspaces just mentioned. 

+ The SpEC code is more cfBcient than the BSSN-OK codes. But SpEC-simulations focus on longer 
and more accurate simulations, resulting in a CPU cost of the same order of magnitude. 
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The vast scven-dimcnsional space of preeessing binaries on circular orbits has 
received so far surprisingly little attention: Campanelli et al [129] compute one 
waveform with mass-ratio q = 1.25, spins of 0.6 and 0.4, lasting about nine orbits. 
A few configurations are used to demonstrate BH-BH mergers in Szilagyi et al. [48], 
but without discussion of gravitational waveforms. Ref. [44] perform several inspiral 
simulations of processing binaries to develop and test eccentricity removal. Sturani 
et al [119,130] perform few-cycle long simulations, exploring a 1-dim line in the 
7-dim parameter space. While these simulations have been used to construct a 
phenomenological waveform- model of preeessing BH-BH [119,130], the short length 
of numerical simulations will severely limit the accuracy of this model. Lousto and 
Zlochower [1] perform 42 simulations to explore black hole kicks with partial spin/orbit 
alignment. Again, the short length of the simulations (about five orbits) makes 
them unsuitable for GW data-analysis purposes. Recently, several papers discuss 
techniques to represent waveforms of preeessing binaries by expanding the waveforms 
in spherical harmonics with respect to a time-dependent frame which is aligned with 
the instantaneous radiation [131-133] (of which [131] presents a preeessing q = 3 BH- 
BH simulation). 

As this short survey illustrates, several groups put their toe into the ocean of 
preeessing binaries, but nobody has seriously braved it yet. The most significant 
current effort to explore preeessing systems is the NR-AR collaboration [122], which 
aims to construct waveform models for preeessing BH-BH systems. While an 
intermediate goal is to revisit the aligned spin case, this collaboration will compute 
several high-quality preeessing waveforms. 

2.4- BH-BH in non-vacuum 

Fully relativistic hydro-dynamics simulations of BH-BH in gaseous environment have 
continued. Bode et al [134] and Bogdanovic et al [135] consider radiatively inefficient 
accretion fiow on merging BH-BH binaries, modeled as a BH-BH embedded in a cloud 
of gas with Gaussian density profile centered on the center of mass. Farris et al [136] 
investigate the BH-BH analog of Bondi-accretion, embedding the BH-BH in ambient 
gas of constant density* . 

Accretion disks are also under continued investigation. Refs. [138-140] study the 
behavior of a circumbinary disk after the BH-BH merger, when the disk responds to 
the remnant black hole with reduced mass (due to energy loss through gravitational 
waves) and with non-zero velocity relative to the center of mass of the accretion disk 
(due to black hole kicks). More recently. Bode et al. [141] and Farris et al. [142] 
investigate BH-BH binaries with a circumbinary disk. 

Two groups are investigating the effects of electro-magnetic fields surrounding 
BH-BH binaries. This work started with solving the vaccum Maxwell equations 
coupled to GR by Palenzuela et al [143, 144] and Moesta et al. [145]. More recently 
a force- free treatment of the electro-magnetic fields was presented in Refs. [146-150]. 
The force- free approximation assumes the presence of a tenuous plasma which shortens 
out any electric field parallel to the magnetic field lines. Such a plasma is produced by 
pair-production and forms the basis for the Blandford-Znajek process [151] to extract 
energy from rotating black holes. Palenzuela et al [146], in particular, demonstrate how 
a merger of a spinning BH-BH can result in a total of three jets: During the inspiral, 

* Zanotti et al [137] studied the standard Bondi accretion onto a single black hole with general 
relativistic radiation- hydrodynamics. 
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one jet associated with each BH; and after the merger, a third jet associated with the 
remnant BH. Moesta et al [150] consider the strength of this beamed emission relative 
to the uniform emission that also accompanies mergers of BH-BH in these cases. 

3. BH NS &: NS NS systems 

We now turn to simulations with at least one Neutron star (NS). These simulations 
have a more varied set of goals than the BH-BH simulations, among them: 

(i) Compute gravitational waveforms to aid GW detectors. This requires simulations 
covering many inspiral orbits at high phase-accuracy, a significant challenge for 
hydro-codes which are typically less accurate than vacuum BH-BH codes. 

(ii) Investigate properties of the merger: In BH-NS binaries, does the Neutron star 
get disrupted? Does an accretion disk form, and how large is it? Is the binary a 
viable progenitor for a short gamma ray burst? 

(iii) Are these mergers a viable source for r-process elements [152]? This requires the 
formation of unbound ejecta which return material into the interstellar medium. 

A variety of physical effects may affect evolution of these binaries and all of these 
effect must be studied numerically: impact of mass-ratio, effect of equation of state, 
effect of magnetic fields, effect of nuclear physics, effect of neutrino transport. For BH- 
NS systems, one must also investigate the effect of the BH spin (both in magnitude 
and direction). Neutron stars are generally assumed to be very slowly spinning, albeit 
some recent work approaches simulations of spinning NS [153, 154]. Several different 
groups investigate BH-NS and NS-NS systems, with each group generally focusing on 
a a subset of these physical effects, as detailed in Sees. 3.1 and 3.2. 

The numerical techniques for BH-NS and NS-NS combine a solver for the Einstein 
equations with appropriate solvers for the matter fields. Most commonly used is the 
BSSNOK evolution system for gravity, combined with high-resolution shock capturing 
techniques for relativistic hydrodynamics, discretized with adaptive-mesh refinement 
finite-differences. The SpEC collaboration combines a spectral treatment of the gravity 
sector (similar to their BH-BH simulations) with a treatment of matter on a separate 
Cartesian finite-difference grid that covers only the region of space in which matter is 
present and that moves along with the Neutron stars. Some recent papers describing 
numerical techniques and individual codes are Refs. [155-163]. Initial data for compact 
object binaries with Neutron stars are described in Refs. [164-169]. 

Detailed investigations into accuracy of BH-NS and NS-NS simulations are given 
in Refs. [170-172]. Broadly speaking, preserving low-phase errors during inspiral 
is difficult and simulations with Neutron stars are less accurate than their BH- 
BH counterparts. No current code is clearly superior; in particular, the advantage 
of SpEC for vacuum simulations does not carry over and SpEC-hydro-simulations 
are of comparable accuracy than other codes. Resolving disk-dynamics is another 
numerically challenging aspect, especially when the disks are long-lived. Finally, MHD 
simulations require a lower cutoff on the matter density. Because gamma-ray-bursts 
involve outflows in regions of low Baryon density, the MHD caveats might become 
important (the IMEX treatment of Palenzuela et al [173] is a recent alternative). 
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3.1. BH-NS 

One key objective of the study of BH-NS binaries lies in finding the region of parameter 
space where the Neutron star is tidahy disrupted. In this case, an accretion disk forms 
and the system is considered a possible candidate for a gamma-ray burst (although 
disruption may not be necessary for a gamma-ray burst, see [174, 175]). Moreover, the 
gravitational wave signature changes dramatically at disruption: Once the Neutron 
star is disrupted and spread into an accretion disk, GW emission is drastically reduced. 

Tidal disruption depends primarily on three parameters. (1) The compaction of 
the Neutron star: more compact stars are harder to disrupt. (2) The mass of the 
black hole: more massive black holes have weaker tidal forces, reducing the tendency 
of disruption. (3) The distance between Neutron star and black hole: For spinning 
black holes, the innermost stable orbit of corotating geodesies moves inward; therefore, 
for BH-NS systems with the BH spin aligned with the orbital angular momentum, the 
orbit of the Neutron star will reach a smaller distance, thus increasing the tendency 
of disruption. 

A significant amount of attention has been devoted on determining the parameter 
space in which accretion disks form and on the properties of the formed accretion disks. 
Let us summarize some recent work: Shibata et al. [176] investigate polytropic Neutron 
stars orbiting non-spinning black holes. For mass-ratios q = Mbh/-^ns = 1-5 and 3, 
disruption occurs, but not so for q = 5. Kyutoku et al. [177] consider BH's with spin 
aligned or antialigned with the orbital angular momentum with spin magnitudes up to 
Xbh = Sbh/M^yi = 0-5- For low mass-ratio q < 3 the Neutron star disrupts. For q — 5, 
disruption occurs only for prograde BH-spin (this work uses a piecewise polytropic 
equation of state, which is fitted to physical equation of states). The remnant disks 
have mass of ~ O.IMq and this work (as the other references mentioned) confirm that 
indeed NS disruption leads to a cutoff in the gravitational wave spectrum. A related 
study by Chawla et al. [178] agrees with the results of [177]: Focusing on mass-ratio 
g = 5, aligned co-rotating BH with spin of xbh = Sbh/M^^ = 0.5, and a polytropic 
equation of state, Chawla et al find that an accretion disk forms, with essentially all 
material gravitationally bound. Ref. [178] also incorporates magnetic fields, and for 
field-strengths up to lO^^G, the effect of the magnetic field was found to be marginal. 

BH-NS binaries with relatively low mass-ratio 9^5 are numerically easier to 
handle and are interesting because these systems form accretion disks easily. However, 
population synthesis suggests that black holes are generally more massive [179]. 
Foucart et al. [171] perform simulations with more massive black holes, AfBH = lOM0 
and mass-ratio q = 7. For a non-spinning black hole, or moderately corotating black 
hole, no accretion disk forms at q — 7. Only at high spins xbh ^ 0.7 does an accretion 
disk develop. These simulations investigate BHs with spins as large as xbh = 0.9, the 
largest to date. 

To close our summary on the effect of BH spin and mass-ratio, we note that 
all simulations mentioned so far have the BH spin parallel to the orbital angular 
momentum. This assumption is removed by Foucart et al. [180] who vary the spin 
direction of the BH. Specifically, the angle 9 between BH spin and orbital angular 
momentum is varied between 0° and 80° (for q = 3, BH spin xbh = 0.5 and a polytropic 
equation of state). Angles 6 > 40° result in a reduction in disk- mass by about a 
factor of 2. For the same mass-ratio, Ref. [180] also investigate aligned BH spins 
with magnitude xbh = 0.9. Consistent with expectations, such a high co-rotating 
BH spin coupled with the low mass-ratio g = 3 results in a very large disk with mass 



Numerical simulations of compact object binaries 13 

approaching ~ 0.4A/ns- 

A second focal point for recent work was the impact of the equation of state 
(EOS) of the Neutron star matter. Duez ct al. [181] investigate polytropic equations 
of state with two different polytropic indices (r = 2 and r = 2.75), and the Shen EOS 
with two treatments of the electron fraction. The binary has mass-ratio q~3 and an 
aligned BH spin xbh = 0.5. Duez et al find that more compact Neutron stars emit 
stronger gravitational waves and result in smaller disk-masses. In particular, tidal tails 
depend on the EOS. Kyotuko et al. [182] investigate piece-wise polytropic equations 
of state [183,184]. They find that less compact Neutron stars tend to disrupt at a 
larger separation, and that disk-mass (and the characteristic GW cutoff frequency) 
correlate with compaction of the Neutron star. 

Recently, the first MHD simulations of BH-NS binaries were performed. Chawla 
et al. [178] consider field-strengths of lO^^G and find no appreciable difference to non- 
magnetic binaries, whereas Etienne et al. [185] report that a magnetic field of lO^'^G 
results in a increased disk-mass (no evidence for collimated outflows were observed, 
although [185] points out that higher resolution simulations would be needed, that 
follow the post-merger dynamics for a longer period of time) . 

Stephens et al. [186] and East et al. [187] consider hyperbolic encounters of a 
Neutron star with a black hole (with relative velocity lOOOkm/sec at large distance). 
Refs. [186, 187] vary the equation of state, impact parameter, and BH spin. By far the 
biggest effect on the results has the impact parameter: Depending on the periastron 
distance, the NS can be disrupted in the first approach, periodic mass-transfer can 
occur, or the NS can pass the BH unharmed. 

3.2. NS-NS 

Binary Neutron stars have seen an equal amount of activity lately as BH-NS binaries. 
The themes are quite similar to BH-NS, with intense efforts of the research groups to 
extend the range of included physical effects. 

Equation of state effects are exhaustively explored by Hotokezaka et al. [188]. 
They simulate six different EOS's (parametrized by piecewise polytropic [183,184]) 
for three different NS-masses each and classify the results into prompt collapse to a 
BH, shortlived hypermassive Neutron star, and long-lived hypermassivc Neutron stars. 
Hotokezaka et al report torus masses of up to 0.1 Mq around the newly formed black 
hole. Sekiguchi et al [189] explore an equation of state with a phase-transition to 
hyperons. This causes substantial differences in dynamics, observable in gravitational 
waves. 

When a NS-NS merger results in a hypermassive Neutron star (as opposed 
to prompt collapse), shocks during the merger will heat the remnant NS to high 
temperatures, and neutrino cooling will become important. Therefore, neutrino 
cooling must be modeled in order to determine reliably the timescale on which the 
hypermassive NS cools and collapses to a BH when the thermal pressure becomes 
insufficient to support the star. Sekiguchi et al. [190] perform the first simulation of 
this process, incorporating neutrino cooling with a leakage scheme, and using a finite- 
temperature Shen EOS. They report the neutrino luminosities for the simulated NS- 
NS mergers, and find indeed that the lifetime of low-mass hypermassive NS depends 
on the neutrino cooling. 

MHD simulations have also been improved. Giacomazzo et al [191] and Rezzolla 
et al [192] focus on very long simulations of the post-merger accretion disk. While 



Numerical simulations of compact object binaries 14 

Rcf. [191] focuses on disk dynamics, Rcf. [192] reports evidence of electromagnetic 
coUimation along the rotation axis of the accretion disk. This is the first claim of 
the "missing link" that connects the post-merger accretion disk with the eventual 
launching of the jets that power gamma ray bursts. Two other research groups 
have published simulations of merging magnetized NS-NS systems [193,194], and 
both these groups have since published improved techniques to handle magnetic 
fields [161,162,173,195]. It will be very interesting to see whether the impressive 
and important results of Ref. [192] can be confirmed by these groups. 

Eccentric NS-NS systems were also studied for the first time. Gold et al. [196] 
find that eccentric mergers result in larger disks, and that periastron passage leads to 
excitation of f-modes in the Neutron stars. 

Finally, recently work has begun to compare NS-NS inspiral simulations with 
post-Newtonian approximations, and to fit analytical waveform models to the 
numerical NS-NS inspiral simulations. Such work is important for gravitational 
wave data-analysis, because complete, phase-accurate waveforms (ranging from early 
post-Newtonian inspiral into merger) promise the most accurate data-analysis for 
gravitational waves from NS-NS. Several NS-NS simulations were performed for these 
purposes all about 10 orbits long. Bcrnuzzi et al. [172] performs a comparison with 
TaylorT4 and finds significant phase-differences (about 1 GW cylce) . which cannot be 
explained by known tidal corrections. Baiotti et al [197, 198] investigate TaylorT4 and 
EOB models. Without free parameters, EOB and TaylorT4 both deviate by about 
one GW cycle from the numerical NS-NS simulation; in contrast, if the EOB model is 
amended with one free fitting parameter (parameterizing the strength of higher order 
tidal effects), the NS-NS simulation can be fitted with an error of only 0.24 radians. 

Read et al. [199] and Lackey et al.[200] analyze large sets of NS-NS simulations 
to determine what information about the equation of state can be extracted from 
future gravitational wave observations. The best-constrained parameter will be the 
compactness of the Neutron star. 

4. Conclusion 

The recent progress in simulations of BH-BH, BH-NS and NS-NS systems is 
spectacular. What are the challenges going forward? 

For BH-BH, the numerical methods are in good shape. Very challenging 
simulations were successfully performed as described in Sec. 2.2, pushing large mass- 
ratio, large spins, and the number of GW-cycles. The open question is whether the 
numerical techniques can handle a combination of these properties, e.g. high spin 
(> 0.9) and high mass-ratio (> 4), and how well accuracy holds up when the number 
of GW cycles in the simulation is increased. Besides these unexplored regions of 
parameter space, the main challenge forward is the systematic exploration of the 
vast parameter space of possible binaries, at sufficient accuracy and length to allow 
gravitational wave detectors to reach optimal sensitivity and optimal accuracy in 
parameter estimation. The computational cost of BH-BH simulations, the possibility 
that the numerical waveforms may have to be much longer than current simulations 
for optimal parameter estimation, and the sheer size of the parameter space will 
necessitate compromises, in length of the simulations or in parameter space coverage 
(or both). The severity of these compromises will only be known after an initial 
exploration of the precessing parameter space and after first attempts to fit analytical 
waveform models. The ease of fitting (currently unknown) will determine how many 
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simulations arc needed. 

Efficiency improvements in the numerical codes will furthermore determine how 
quickly the parameter space can be sampled and how quickly on-demand simulations 
can be performed (e.g. in response to gravitational wave detectors observing a 
tentative event). Such efficiency improvements may come from novel computer 
algorithms like implicit time-stepping [109,110] or from utilizing novel computing 
paradigms like graphical processing units. 

For BH-NS and NS-NS systems, the current frontier is exploration of all relevant 
physical effects. Given the complexity of the simulations and the varied micro-physics, 
it is imperative that different groups perform similar simulations to cross-check. Once 
qualitative features are explored, the field will turn towards quantitative sampling 
of the parameter space, with systematic and careful calculation of gravitational 
waveforms. 

Numerical and implementation issues of compact object binaries seem a problem 
of the past, as attested by the stunning progress in numerical relativity. The deeper 
insight into formulation of the equations, and experience of what works and what 
does not makes it now possible to consider entirely different numerical algorithms, 
like discontinuous Galerkin methods [201,202], moving Voronoi meshes [203,204], or 
novel time-stepping techniques [109]. 
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